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Abstract 

A new approach is presented to evaluate multi-loop integrals, which appear in the calculation of cross-sections in 
high-energy physics. It relies on a fully numerical method and is applicable to a wide class of integrals with various 
mass configurations. As an example, the computation of two-loop planar and non-planar box diagrams is shown. The 
results are confirmed by comparisons with other techniques, including the reduction method, and by a consistency 
check using the dispersion relation. 
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Si 

q. 1. Introduction 

In the study of high-energy reactions observed at current and future accelerators, such as LHC and ILC, precise 
theoretical predictions of cross-sections including higher order corrections are required. This is due to the fact that 
the lowest order approximation in perturbative calculations of quantum field theory is not sufficiently accurate to be 
compared to the experimental data. One has to take into account the contributions from higher order terms as well. 
In order to include these corrections in the Standard model or beyond, it is indispensable to handle the evaluation of 
loop integrals. 

At the one-loop level it is known that analytic solutions exist for any type of diagram, and the results are expressed 
in terms of known functions, such as logarithms and Spence functions (see, for example [1]). Using these analytic 
results several automatic computation systems [2, 3, 4, 5, 6, 7, 8, 9, 10] have been proposed. In order to estimate 
cross sections we need automatic computation systems because we may have to deal with a large number of relevant 
Feynman diagrams for a given process. 

However, the extension of the system to include higher order corrections is not an easy task, because analytic 
integration is generally impossible for higher loop diagrams, especially for diagrams which depend on more general 
mass configurations. Analytic results are only known for a limited class of two-loop diagrams. Therefore we have 
to rely on numerical evaluations. We need to establish efficient methods that can be incorporated into automatic 
computation systems of cross-sections. For a number of years we have gained experience evaluating one-loop integrals 
numerically, where the results can be compared with known analytic answers. We succeeded in calculating vertex, 
box and pentagon diagrams with arbitrary masses [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. We also computed 
two-loop self-energy and vertex diagrams. Further related work can be found in [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 
33, 34, 35, 36, 37, 38, 39]. 

In our method we start from the Feynman parameter representation of loop integrals. We employ a fully numerical 
integration procedure combined with numerical extrapolation. The purpose of this paper is to describe the method in 
detail and to show results for more complicated loop integrals, corresponding to two-loop box diagrams with massive 
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particles. For simplicity we deal with scalar loop integrals throughout this paper, ignoring all spin complications that 
are irrelevant to the essential discussion of the numerical approach. 

The most general form of the scalar integral for a diagram with L loops and N internal lines is given by 



N 



j — l x ' r— 1 

where lj is the j-th loop momentum in the n-dimensional space-time, and 

D r = ql - ml + ie (2) 

is the inverse of the r-th Feynman propagator, where e denotes an infinitesimal quantity, m r is the mass of the r-th 
particle, and the momentum q r flowing on the r-th internal line is given by a sum of loop and external momenta. We 
make use of the Feynman identity, 

Carrying out the loop momentum integrations delivers 

"(±P-("-t)"'. 

where 

fl N ffN-n(L+l)/2 

I = { -^ N l T\ dx ^ 1 -T,^ {D _ ieC) N- nL/ , - (5) 

The function D is a polynomial in the Feynman parameters {xi}. D further involves physical variables such as the 
external momenta and particle masses. The function C is also a polynomial in the {xi}. Both functions are determined 
by the topology of the Feynman diagram. Details of their construction are summarized in Appendix A. 

In the two-loop box diagrams, D depends on two kinematical variables s and t, where s is the square of the total 
energy of the colliding particle system, and t < is the squared momentum transfer between the initial and the final 
particles. For the infrared divergent integrals, we have two prescriptions. One is to introduce a small fictitious mass A 
for the massless particles and the other is the dimensional regularization technique. In the former we can set n = 4 in 
Eq.(5) and the procedure is straightforward once the value A is fixed [14, 16]. For the latter we put n = 4 + 25 = n{5) 
and use a double extrapolation technique for both e and 8 (from n(6)) in Eq.(5) [20, 21]. Here we estimate the integral 
for a fixed value of S using the extrapolation with respect to e. Repeating this for a series of S values, we can estimate 
the pole residue of 1/6 and the finite part of the integral numerically. 

We briefly describe the general properties of the integral I and give some terminology. Depending on the value of 
s, the function D in the denominator may vanish in the integration domain. In this case, the infinitesimal parameter e 
prevents I from diverging. Then I exhibits an imaginary part even if all the physical parameters s, t and the masses 
are real. This region of s is called the physical region, where s exceeds the threshold energy, so that the reaction takes 
place. On the other hand, in the unphysical region, s is lower than the threshold. This is the region of s where we 
can put e = and the integral is real for real s. Thus the integral I can be regarded as an analytic function in the 
complex s-plane with cuts along the real s-axis, starting at branch points which are determined by physics conditions. 
However, as we shall see below in Section 3, we treat e not as infinitesimal but as a finite number in the numerical 
procedure for calculating 3?(7) and I = 5R(7) + i^s(I), in the physical region. 

This paper is outlined as follows. In Section 2 we construct the integrands for two-loop box diagrams (L = 2 
and N = 7), and present suitable variable transformations. We explain the details of our techniques in Section 3; 
and the results of the computations are shown in Section 4. Section 5 is devoted to a discussion on how to assess the 
correctness of the obtained results. Section 6 gives conclusions and future directions for this work. 
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2. Two-loop box integrals 

The topology of the two-loop box diagram is depicted in Fig. 1. We call Fig. 1(a) the planar diagram and Fig. 1(b) 
the non-planar diagram, respectively. The loop integral in the Feynman parameters (x\, ■ • - , x 7 ) is of the form 



,1 7 

I = — I dx\ dx2 dx 3 dxi dx^ dxe, dx 7 8(1 — X() 
Jo 



C 



(D - ieCf 



Here, D and C are polynomials of Feynman parameters. Their derivations are given in Appendix A. 
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Figure 1 : (a) Two-loop planar box diagram (b) Two-loop non-planar box diagram 



The external momenta pi,P2, P3 and are defined to flow inward, satisfying pi+P2+P3+Pi = 0. The kinematical 
variables s and t are given by 

s = (Pi +P2) 2 = (p3 +Pi) 2 , 
t = (P1+P3) 2 = (P2+Pi) 2 - 
For later notational convenience we introduce a third kinematical variable u by 

u = {pi +Pi) 2 = (P2 +P3) 2 - 
The variables s, t and u are not independent, as 

s + t + u=p\+p\+pl+ p\. 

In the following we derive the explicit formulae of the functions D and C. We also show examples of variable 
transformations, which allow eliminating a common factor in the numerator and denominator. Furthermore, the re- 
sulting form of the integral will be suited for an application of the reduction formalism given in Section 5.1. Followed 
by a Monte Carlo integration we will use the latter for the purpose of comparing of its numerical results with those by 
DCM. 

2.1. Explicit formulae of the functions D and C for the planar diagram 

The functions D and C in Eq. (6), corresponding to Fig. 1(a), are given by 

D = C^xprn} (7) 

- {s(xiX 2 [Xi + X 5 + Xe, + Xj) + X 5 X 6 (xi + X 2 + £3 + X4) + XiXiXe, + X2X4X5) 

+ tX3X4Xj 

+ Pi(x3{xiX4 + XiX 5 + XxXe, + + X4X5)) 

+ pl(x 3 (x 2 X4 + X 2 X 5 + X 2 X 6 + X 2 X 7 + X^Xq)) 

+ pl(x 7 (xiX4 + XlX 5 + X 2 X 5 + X 3 X 5 + 

+ pl(x 7 (xix e + X 2 Xi + X 2 Xe, + X 3 Xq + X^Xe))}, 
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and 



C = (xi + x 2 + x 3 + X4)(x4 + x 5 + x e + x 7 ) - x\. (8) 

Using the transformation (xi, X2, x 3 , X4, x$, x§, £7) — > (pi, p 2 , p 3 , u\, u 2 , u 3 , U4), defined by x\ — p\U\, x 2 — 
piu 2 , x 3 = pi(l - u\ — u 2 ), x 4 = p 3 , x 5 = p 2 u 3 , x 6 = p 2 U4 and x 7 = p 2 (l — u 3 — u 4 ), we obtain 

J dxi ■ ■ ■ dx-j <5(1 — ^ Xj) = J dpi dp 2 dp 3 6(1 — ^ pj) p\p\ J dui du 2 du 3 du 4 , 



with 



Pi + Pi + P3 = 1, < u\ + u 2 < 1, < u 3 + u 4 < 1, 



and with Jacobian p\ p\. Changing the variables by pi = p£, p 2 = p(l — £) and p 3 = 1 — p gives 

J dp x d P2 d P3 s(i - j2 Pi) p\p\ ■ ■ ■ = I d P J o <% P 5 ea - o 2 • • • • 

After these transformations, D and C contain a common factor p and we set V = D/ p and C = C/p. The integral 
becomes 



Iplanar 



dp dS, j dui 
Jo 
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du 2 / du 3 
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du 4 
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(V - ieCf 



p 6 e(i-zr, 



where V is a quadratic in u = (m, u 2 , u 3} w 4 ) T , 
V = u T Au + B T u + c, 

and 

C = rtil - + 1 - p. 



(9) 



(10) 



(11) 



The 4x4 matrix A is symmetric and depends on the internal masses mi (1 < £ < 7), and on the kinematical 
variables, s and t. In this paper we assume p 2 = p 2 = p 2 = p 2 = m 2 for both diagrams. When mi = m 2 = to 5 = 

m & — m and m 3 = to 4 = m 7 = M (to 7^ M), we have 

A _( P e(i-pOAi P t(i - P )(i - oa 2 \ 

U(! - - 0^2 P(l - 2 (1 - P(l - W ' 



where the 2 x 2 matrices Ai and A 2 are 



— m 2 s/2 — w 
s/2 - to 2 



t/2-m 2 s/2 + f/2- 

s/2 + t/2-m 2 t/2-m 2 



The vector B is given by 



B 



/ -tp£(l-p)(l-0+M 2 ptC \ 

-tpi(i-p)(i-i) + M 2 P zc 

-t P ai-p)(i-o+M 2 p(i-oc ' 

\-tpt(i-p)(i-0 + M 2 p(i-(;)cJ 

and the scalar c is 

c = tpt(l-p)(l-0-M 2 C. 
The quadratic form will further be used in Section 5.1 for a comparison of DCM with a reduction method. 
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2.2. Explicit formulae of the functions D and C for the non-planar diagram 
The functions D and C in Eq. (6), corresponding to Fig. 1(b), are 

D = C'y*^xem% 

- {s(xiX 2 X i + X\X 2 X 5 + XlX 2 X 6 + XlX 2 X 7 + XlX 5 X 6 + X 2 X 4 X 7 - £3X4X6) 

+ i(x 3 (-X 4 X 6 + X5X7)) 

+ P 1 (x 3 (X\X4 + X1X5 + XlX 6 + X1X7 + X 4 X 6 + X4X7)) 

+ pl(x 3 (x 2 x i + x 2 x 5 + x 2 x 6 + X 2 X 7 + X 4 X 6 + X 5 X 6 )) 

+ pl(xiX 4 X 5 + X1X5X7 + X2X4X5 + X 2 X 4 X 6 + X3X4X5 + X3X4X6 + X4X5XQ + X4X5X7) 

+ _P4(xiX 4 X 6 + XlX 6 X 7 + X2X5X7 + X 2 X 6 X 7 + X 3 X 4 X 6 + X 3 X 6 X 7 + X 4 X 6 X 7 + X 5 X 6 X 7 )} 

and 

C = {x\ + x 2 + x 3 + x 4 + x 5 )(xi + x 2 + x 3 + x 6 + x 7 ) - (xi + x 2 + x 3 ) 2 . 

The transformation (xi, x 2 , x 3 , X4, x 5 , x 6 , X7) — > (pi, p 2 , p 3 , u\, u 2 , u 3 , U4), defined by x\ — piU\, x 2 
x 3 = pi(l - m - u 2 ), x 4 = p 2 u 3 , x 5 = ,02(1 - u 3 ), x 6 = p 3 M 4 and x 7 = ,03(1 - u 4 ), yields 

J dxi ■ ■ ■ dx 7 (5(1 — ^ Xj) = J dpi dp 2 dp 3 S(l — ^ pf) p\p 2 p 3 J dux du 2 du 3 du 4 , 

with 

Pi + Pi + P3 = 1, < u\ + u 2 < 1, < u 3 < 1, < U4, < 1, 
and the Jacobian is p\p 2 p 3 . The change of variables p\ = 1 — p, p 2 = p£ and ,03 = p(l — £) gives 

J d Pl dp 2 d P3 6(i-j2pj)plp2p3--- = I dpJ o ^P 3 (i-p) 2 e(i-o--- , 

with 

0<p<l, < £ < 1. 

Similar to the case of the planar diagram, D and C contain a common factor p. Putting V = D / p and C 
delivers the final form of the integral 

Inon-planar = ~ / dp d£ du x \ du 2 \ du 3 \ du 4 — — — T P(l - pf^ ~ , 

Jo Jo Jo Jo Jo Jo \ U ~ leL ) 

where I? is a quadratic in u = (u\, u 2 ,u 3l u 4 ) T 1 given by 
V = u T Au + B T u + c, 

and 

c = P i{i - + 1 - p. 

With the mass assignment as mi = m 2 = m 4 = m 6 = m and m 3 = m 5 = ra-j = M(m ^ M) we have 

A _f (1-P) 2 A! tf(l-p)(l-QA 2 

^(l-p)(l-0^2 A 3 

where the 2 x 2 matrices A\, A 2 and A 3 are 



s/2-m 2 \ _( t/2-m 2 s/2 + t/2-m 2 

Al -\ S /2-m 2 -m 2 J' A2 ~ \s/2 + t/2 - m 2 t/2-m 2 
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A 3 



( 



(- S /2 + m>(l-p)£(l-£) 



(- S /2 + m 2 )p(l-p)ai-0 
-m 2 p(l-C) 2 (l-p(l-0) 



') 



The vector B is given by 



B = 



/ -tpC(l-p)(l-0+M 2 (l-p)C\ 
-tpt(l-p)(l-0+M 2 (l-p)C 

-tpai-p)(i-a)+M 2 ptic 

\-tpZ(l-p)(l-Z) + M 2 p(l-Z)Cj 



and the scalar c is 



c = tp(l-p)t(l-0-M 2 C 



Note the similarity between the Eqs. (9) and (14) of the planar and the non-planar integral, respectively, obtained 
via the transformations in Sections 2. 1 and 2.2. The transformed integrand functions (of both 6-dimensional integrals) 
involve a function V in the denominator, which is a quadratic in the variables u\, u 2 , u 3 and w 4 . This form of the 
integrand will further lend itself to the reduction method of Section 5.1 (which will be used for a comparison of the 
numerical results). As an aside, the form of the 4-dimensional integral in u\, u 2 , u$, U4 also resembles that of the 
one-loop pentagon integral, e.g., in [19]. 

3. Numerical techniques 

We introduce the Direct Computation Method (DCM), based on a combination of numerical integration, and 
extrapolation on a sequence of integrals. DCM comprises the following three steps: 

1 . Let e in Eq. (5) be a finite value determined by a (scaled) geometric sequence 



for a constant e and base < l/A c < 1. 

2. Evaluate the multi-dimensional integral / of Eq. (5) numerically. In view of the finite ei we obtain a finite value 
for the integral corresponding to each I. Thus a sequence of I(ei), I = 0, 1, 2, • • • is generated. 

3. Extrapolate the sequence I(ei) to the limit as e; — > with the purpose of calculating I as lim e ^ ^(e)- 

If D does not vanish within the integration region, we can ignore e and no extrapolation is needed as / = 1(e) | e=0 . 

For multi-dimensional integration we make use of the DQAGE routine in the QUADPACK [40] package. DQAGE 
uses a variant of Gaussian quadrature, where the sampling points are given by a Gauss-Kronrod rule pair in each 
subinterval. The Gauss rule with v points has polynomial degree of accuracy d v = 2v — 1; i.e., it is exact for 
polynomials of degree d = 0, ■ • ■ ,d u and not for all polynomials of degree d v + 1. The corresponding Kronrod rule 
re-uses the abscissas of the Gauss rule and adds v + 1 points interlacing with those of the Gauss rule. The Kronrod 
rule with 2^ + 1 points has polynomial degree 3^ + 1 if this number is odd (for v even), and otherwise 3^ + 2 (for v 
odd). 

On input for DQAGE, the user selects one of six Gauss-Kronrod pairs, with 15, 21, 31, 41, 51 or 61 points, via 
the input parameter key = 1, 2, 3, 4, 5, or 6, respectively. The rule pair produces the Kronrod rule value as the 
integral approximation, together with an estimate of the absolute error on each subinterval (which is based on the 
difference between the Gauss and the Kronrod result on the subinterval). This allows the selection of that subinterval 
with the largest estimated error, as the next interval to be subdivided in successive steps of the adaptive partitioning 
strategy of DQAGE. The user imposes a bound on the number of subdivisions via the input parameter limit. As a 
result of the adaptive partitioning, the algorithm subdivides intensively around singularities, so that hot spots emerge 
where singularities or other irregular integrand behavior occur within the integration interval. For multi-dimensional 
integration we apply DQAGE in a repeated (iterated) quadrature for successive coordinate directions [41]. 

In DCM, the accuracy of the result depends on that of the calculated sequence of integrals /(e;), I = 0, 1, ■ • ■ . 
Since the integration error affects the accuracy of the extrapolation, we want to compute the integrals I(ei) to at least 
an order of magnitude more accuracy than that expected for the final result. On the other hand, the CPU time is 



e = e; = e /(A c ) 1 , 1 = 0, 1, • • • , 



(17) 
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directly related to the accuracy requirement. Thus, adequate values need to be specified for the input parameters key 
and limit of the DQAGE routine in each dimension, to control the overall work and the termination of the algorithm. 
For the computation of two-loop box integrals, we find that key = 1 or 2 and limit = 10 ^ 30 are suitable values. 

We use Wynn's e-algorithm[42, 43] for the extrapolation, which works efficiently under fairly general conditions, 
even for very slowly convergent sequences or series. The e algorithm is applied to the sequence I(ei), I = 0, 1, • • • 
obtained by multi-dimensional integration. We define the table elements a(l, k) of the extrapolation table with initial 
values 

a(/,-l) = 0, a(l,0) = I(ei), f = 0, 1, - - - . (18) 
The element a(l, k + 1) is obtained from a(l, k) and a(l,k — 1) by the following recurrence relation: 

<^ + l)=a { l + l,k-l) + a{l + J ) _ a{hky 1 = 0,1,-... (19) 

Whilst the a(l, fc)'s with odd k are meant to store temporary numbers, the a(l, k)'s with even k give extrapolated 
estimates. 

We use the e-algorithm code from the QUADPACK [40] package. With each new I{e{), a new lower diagonal can 
be added to the extrapolation table. At each iteration only the last two lower diagonals need to be stored for this 
computation. Along with each new table element a{m, n) where n is even, an error estimate is calculated based on 
differences with its neighboring elements. In a converging table, the even-numbered columns as well as the diagonals 
converge to the limit lim e ^ I( e ) = ^ (barring roundoff). Among the even-column indexed table elements along 
newly computed lower diagonal, the e-algorithm code selects the "best" a(m, n) (with the least error estimate). The 
CPU time for the extrapolation is negligible compared to that of the integration. 

We further have to use some heuristics for the computation of the extrapolated sequence. The acceleration constant 
A c in (17) can usually be set to 2. In cases where the integration is very difficult for decreasing e; , a smaller value of 
A c is used, e.g., A c = 1.3 or 1.2, yielding a sequence of e;, I = 0, 1, • ■ • which decreases more slowly. To determine 
the initial value of the geometric sequence, we assign e depending on the squared mass appearing in function D. We 
parametrize e in the form e = AJ, where the parameter 7 can be adjusted. The choice of these parameters influence 
the accuracy of the result. 

For the two-loop box integral computations reported in the next section, we found A c = 1.2 and 7 around 40 to 
be adequate values. The accuracy achieved is restricted by the actual CPU time needed. If the computation time is 
excessive, we have to accept less accurate results. This happens, for example, when s is much greater than 10m 2 . 
Thus the accuracy is different from point to point in the plots shown below. All the computations are done in double 
precision arithmetic. 



4. Numerical results 

According to the prescription of DCM in the previous section, we evaluate both I p i an ar and I n on-planar given by 
Eqs. (9) and (14), respectively. In both cases the kinematical variable s is varied but t is fixed at t = -10000GeV 2 
throughout the computations. We introduce the dimensionless variable 

fs = ^- (20) 
For the mass parameters we set m — 50 GeV and M = 90 GeV. 



4.1. Planar diagram 

In previous work [19] we presented results of the real part integral in the physical region, 4.5 < f s < 25.0. Here 
we evaluate the integral in the region 0.0 < f s < 25.0 for the real part and the imaginary part. The results are depicted 
in Fig. 2 where the data points represent the integral values and the lines merely connect the points as a guide for the 
eyes. The s-channel threshold starts at f s = 4.0 corresponding to s = 4m 2 . For example, we set e = 1-2 45 , key = 1 
and limit = 10 in all dimensions for the real part at f s = 10.0 and it took 8.5 days to obtain the result with enough 
accuracy as 0.01% using a system with Intel Xeon CPU E5430 @ 2.66GHz. 
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4.2. Non-planar diagram 

Fig. 3 shows the results for —20.0 < f s < 20.0. Different from the planar case, it is known that I n on-planar has 
two cuts; one starts from the normal s-channel threshold, s = Am 2 , and the other from s = — t — M 2 — AmM to 
s = — oo. The latter corresponds to the w-channel threshold at u = (M + 2m) 2 . These correspond to f s = 4.0 and 
f s = —6.44, respectively. In Fig. 3 we also show some results of the imaginary part in the range —100.0 < f s < 
—20.0. In this region the imaginary part is small but its contribution is not negligible when it is put in the dispersion 
integral (32) of Section 5.2. 

In the physical region the computation time tends to be longer for larger f s . This applies to both the planar and the 
non-planar diagram. For example the time required to obtain the real part of the non-planar box integral with enough 
accuracy as 0.003% at f s = 10.0 (with e = 1.2 40 , key = 2 in all dimensions, and limit = 10, 20, 20, 10, 10, lOin 
consecutive dimensions) is about a week using a system with Intel Xeon CPU X5365 @ 3.16GHz. For much greater 
f s it may become more difficult to get an answer in a practical time. However, this computation time is measured 
using a single CPU. It can potentially be shortened by applying parallel computing techniques on (possibly distributed) 
multi-core processors [17, 44, 45]. 



5. Validation of the results 

After obtaining answers by the numerical computation, the most important issue is how to confirm that the results 
are correct and reliable. It would be most desirable to have answers available by independent methods. For example, 
we computed the result 0.10364072096 ± (0.315 x 10~ 7 ) for the planar integral, with s = t = 1 GeV 2 and m\ = 1 
GeV 2 (1 < i < 7), v\ = 1 GeV 2 (1 < i < 4). We were able to compare this to the value 0.1036407209893 
evaluated by the program SYS [46] and found good agreement. With the same values of the kinematical variables we 
obtained the non- planar integral as 0.08535139±(0.105 x 10~ 7 ), but no result was available by SYS. This comparison 
demonstrates that the expressions of the functions C and D for the planar diagram are correct and that DQAGE works 
as expected. In these examples, D does not vanish within the integration region; thus we do not need extrapolation. 
The CPU time required for both the planar and the non-planar diagram is less than 2 min. using a Xeon CPU X5365 
@ 3.16GHz. 
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Below we outline an integration method based on reduction formulas, and explain how to use it for a consistency 
check. 



5.1. Comparison with the reduction method 

Consider a quadratic form V in N variables, u = (m, ■ ■ ■ , ujv), 

V = u T Au + B T u + c, (21) 

where A is an iV-dimensional symmetric matrix A = (Aij), B is an TV-dimensional vector B T = (Bi, ■ ■ ■ , Bfq), 
with constant coefficients. Here c contains not only a real part but also — ieC with an infinitesimal e. Assuming A is 
invertible, we define the vector 

X T = 2u T + B T A-\ (22) 

Then we have 

X T (V£>) =4(V-c)+B T A- 1 B = 4V + A N , (23) 
with V T = (d/dui, • • • , d/dujsi)- Here A^ is defined as 

A N = B T A" 1 B - 4c. (24) 

We divide Eq. (23) by V a+ 1 where a > is an arbitrary number. Using the relation 

X \ 2N X T VP 

- (25) 



we obtain the following reduction formula 

A„ _-4 + 2W/o _l vI /X a 



pa a V 25 

It should be noted that the power of the denominator in the right-hand side is decreased by one, compared to the 
left-hand side, that is, the singular behavior is softened. When a — we find 

^ = -4-2iVlog£> + V T (Xlog2?). (27) 

When a polynomial in u, /(u) ^ 1, occurs in the numerator of the left-hand side, the formulas are generalized to 

^ = -4/ + Vr ( /X)/a_l /gey 

and 

/Aw 



= -4/ - V 1 (/X) log 2? + V 1 (/X log V) . (29) 



We apply the formula to the functions V given in Eqs. (10) and (15), which are quadratics in u\, ■ ■ ■ ,114. Since we 
have TV = 2a (N = 4, a = 2) in both cases, we find a simpler formula 

Al " W*Y (30) 



£> 3 2 V^ 2 7 
with A 4 = A 4 (/9, £). By integrating we have 
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Figure 4: Numerical results of ?R(I n on-planar) in units of 10 12 GeV for 0.0 < f s < 20.0. Values calculated by DCM and those by RM 
are shown by circles and bullets with error bars, respectively. 

In Eqs.(9) and (14), the above expression is integrated over p and £. Here A4 = can occur in p-£ space, and is 
regularized numerically by setting the integrand to zero in the vicinity of this anomaly. 

The right-hand side is immediately integrated once. Applying the reduction repeatedly to the form in the right- 
hand side, we see that the original integral is finally replaced by a sum of integrals of functions involving logarithms. 
Thus the severity of the integrand singularity is reduced, which allows performing the integration even with Monte 
Carlo routines. Note that this procedure generally creates lengthy expressions. The imaginary part results from the 
logarithms; let R be a positive number and let z — —R ± ie, then log z — \og(—R ± ie) = log R ± in. We refer to this 
integration method as the Reduction Method (RM). We computed the two-loop box integrals by using BASES [47]. 
The real part of the planar diagram integral in the physical region, shown in [19], is in good agreement with the results 
by DCM. On the other hand, for the imaginary part, the Monte Carlo integration failed to convergence satisfactorily. 

In Fig. 4 we show the real part of the non-planar case obtained by the reduction formulas. Agreement with the 
results by DCM is poor around the threshold f s = 4 in view of poor convergence of the integration by RM. This may 
be caused by the numerical regularization in the vicinity of A4 = as mentioned above. 

5.2. Consistency check using dispersion relation 

The dispersion relation provides a good tool for a consistency check. Based on the observation that I(s) can be 
regarded as an analytic function in the complex s-plane, the real part and the imaginary part, $l(I(s)) and 3(/(s)), of 
J(s) satisfy the dispersion relation, 



where P denotes principal value integral. Recall that I(s) is real in some region of s and accordingly 3(/(s)) vanishes 
there. This integral relation, which is the consequence of the analyticity of I(s), implies the real part can be estimated 
from the imaginary part. DCM computes the real and the imaginary part independently, as they are given by separate 
integrals. However, the dispersion relation indicates that both parts are not independent. They should be consistent 
with the relation of Eq. (32). 




(32) 
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Figure 5: Numerical results of $t(I n on— planar) m units of 10~ 12 GeV - 6 for 0.0 < f 3 < 20.0. Values calculated by DCM and those by the 
dispersion relation are shown by circles and squares, respectively. 

In order to show how this relation works we consider the two-loop non-planar box as an example. In this case the 
integral can be written as 

where so = 4m 2 and s' Q = — t — M 2 — AmM are the threshold in the s-channel and in the u-channel, respectively, as 
described in Section 4.2. 

For the principal value integral computation we used the trapezoidal rule, assuming that 3(/(s)) = far away 
from the origin, for f s < —100.0 and f s > 50.0. Values of 5R(7(s)) resulting from this computation are plotted in 
Fig. 5 for 0.0 < f s < 10.0. The results show good agreement with those by DCM. Thus the relation of Eq. (32) 
enables a consistency check for the answers produced by DCM. 

6. Conclusions 

In this paper we calculated the scalar integrals of two-loop planar and non-planar box diagrams involving massive 
particles. We introduced the Direct Computation Method (DCM) for the evaluation. The novel idea in DCM is that the 
e value in the propagators is treated numerically as a finite number, not as an infinitesimal value. In view of the finite e. 
the integrand of the loop integral is no longer singular. The integration can be carried out numerically for both the real 
part and the imaginary part. Consecutive integrations, for each e/ , produce a sequence of integrals I(ei ) , I = 0, 1 , • • • , 
which are supplied to the extrapolation procedure. A numerical answer for the loop integral results in the limit as ei 
tends to 0. 

Since DCM does not impose restrictions on the values of mass parameters, the method is valid when masses are 
complex [17]. For this case we can put e = where no extrapolation is needed in the same manner as in the non- 
physical region. This flexibility is remarkable and useful for the calculation of cross-sections where decaying particles 
are involved. 

In order to check our evaluations we compared the results with those obtained by other methods [46], including 
Reduction Method (RM). Comparisons of the results have shown satisfactory agreement. The examination of the 
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dispersion relation has lead to a consistency check between the real and the imaginary part of the integral. Thus we 
have established various ways to confirm the results. 

Some issues linked with the implementation of parameters, e.g., the choice of e; values, need to be solved heuris- 
tically. Furthermore, in some regions of the kinematical variables, DCM requires very long CPU times to obtain 
reasonable accuracy. It may be possible to tackle the CPU time problem by utilizing recent developments in computer 
resources and parallel computing technologies [17, 44, 45]. Throughout this paper we use double precision arithmetic, 
but quadruple or extended precision may be needed for some mass configurations, including a small fictitious mass A 
to regularize infrared divergent integrals [16]. This can be incorporated in dedicated program packages [48, 49] . 

For a specified high-energy reaction, all the necessary two-loop diagrams can be generated automatically using 
the GRACE system [4]. The next stage, which involves the automatic generation of amplitudes (i.e., the integrands of 
loop-integrals), would be manageable in view of the experience we gained in handling tree and one-loop processes [4]. 
Thus the only component which needs further development for the construction of an automatic computation system 
for two-loop reactions is a robust loop integral evaluation system. 

Concerning the further development of DCM we need to test integrals for various mass configurations different 
from those in this paper, particularly, infrared divergent integrals by the prescription using a fictitious mass. We 
also need to examine loop integrals with a non-trivial numerator, and explore a systematic treatment of ultra-violet 
divergence. From a technical point of view, reducing CPU time and automatic tuning of the integration parameters 
should be included. After completion of these studies, we expect that DCM will play an important role in constructing 
automatic computation systems for higher-order corrections. 
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Appendix A. Construction of the functions D and C 

For a given diagram, the explicit form of the functions D and C is determined by the following steps [50, 5 1 , 52] . 
• Step 1. 

N 

1. Assign the parameter Xi to the i-th internal line. The parameters {x^ satisfy Xi = 1. 

i=l 

2. Define L topologically independent loops and label them as a = 1, • ■ • , L. The loop momentum l a flows 
through the a-th loop in its own direction. 

3. External momenta pj, j = 1, • • • , K are presumed to enter the diagram inward. Here K is the number 
of external lines. We let pj flow through the diagram while respecting the momentum conservation at 
each vertex. A simple example is where each pj, j — 1, • • • , K — 1 flows through the diagram along a 
continuous path, to reach the vertex where px enters. In this case the momentum conservation is trivial 

4. Each internal line has its direction and the momentum ki for the i-th internal line is defined along this 
direction. It can be expressed by a linear combination of the l a and pj as 

L K 

a=l j=l 

where 

{1 l a flows along the i-th internal line parallel to its direction 
— 1 l a flows along the i-th internal line anti-parallel to its direction 
(otherwise) 

and <t] can be defined in a similar manner fox pj. We define p ex t,i — Y^j=i ^jPj f° r the i-th internal line. 
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5. It should be noted that, even though the choice of the L-loops, the selection of loop-momenta l a , and the 
flow of external momenta are not unique, the final result is the same for any choice. 

• Step 2. 

We construct an L x L symmetric matrix A, an L-vector B and a scalar c. 

N N N 

A a b = ^2<T l a <T l b Xi, B a = ^2<jlx l p ext ,i, c = ^Xi(pl xt i - m\). 

i—1 i—l i—1 

• Step 3. 

The functions C and D are obtained by 

C = det(A), and D = -dot ^ £ T ^ 
D is a homogeneous polynomial of degree L + l, and C is a homogeneous polynomial of degree Lmxi. 

Appendix B. Two-loop box diagrams 

Following these prescriptions, D and C (Eqs. (7), (8)) for the two-loop planar diagram (Fig. 6(a)) are obtained 
from 

An = xi + x 2 + x 3 + x 4 , A12 = A 2 i = —x 4 , A 2 2 = x 4 + x 5 + x e + x-j, 
Bi = xipi - x 2 p 2l B 2 = x 5 pi - x 6 p 2 + x 7 (pi + p 3 ), 
c = x-i(pl - m\) + x 5 (pj - m 2 5 ) + x 2 {p\ - m\) + x 6 (pl - ml) 
+x 7 ((pi +P3) 2 - m%) + x 3 (-m 2 3 ) + x 4 (-ml), 

and those (Eqs. (12), (13)) for non-planar diagram (Fig. 6(b)) from 

An = xi+x 2 + X3 + x 4 + x 5 , A12 = A 2 i = xi + x 2 + x 3 , 
A22 = Xi+X 2 + x 3 + x 6 + XT, 

Bi =-(xi+ x 4 )p 3 - x 3 (p 1 + p 3 ) + x 2 p4, B 2 = (x 2 + x e )p 4 - x 3 (px + p 3 ) - xip 3 , 
c = x x {p\ - m\) + x 4 (pl - mf) + x 2 {p\ - m§) + x 6 (p 2 4 - m%) 
+x 3 {{pi +p 3 ) 2 - ml) +x 5 (-ml) +x 7 (-m%). 



14 



Pi 



P-2 




x 3} m 3 



x 2 ,m 2 



x e ,m e 



Pi P2 



x 2 ,m 2 



x 6 ,m 6 



(a) 



(b) 



Figure B.6: The quantities in Appendix B are obtained from the configuration shown in the figure for (a) two-loop planar box and (b) two-loop 
non-planar box. The arrow on each internal line defines its direction(Step 1.4). Red lines and blue lines show the flow of loop momenta(2 a , Step 
1.2) and that of external momenta(pj, Step 1.3), respectively. 
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